| Descargar notebook.

4. Estimación por Kriging (UG 4)

Nuestro objetivo es estimar el contenido del commodity de interés utilizando solo aquellas muestras pertenecientes a una misma población (UG) en el volumen espacial.

De momento, no sera posible limitar la estimación dentro del limite definido para cada una de las Unidades Geológicas en la sección anterior, y esto quedara para la sección siguiente como post-proceso una vez finalizada la estimación en cada UG.

Se entrega aca el flujo de trabajo a modo de ejemplo sólo para una unidad geológica (UG 4). Usted y su grupo debe replicar el procedimiento para el resto de las unidades.

4.1. Cargar las bibliotecas requeridas

El siguiente código carga las bibliotecas requeridas.

import numpy as np                                            # ndarrays for gridded data
import pandas as pd                                           # DataFrames for tabular data

import geostatspy.GSLIB as GSLIB                              # GSLIB utilities, visualization and wrapper
import geostatspy.geostats as geostats                        # GSLIB methods convert to Python
import geostatspy
print('GeostatsPy version: ' + str(geostatspy.__version__))

import os                                                     # set working directory, run executables

from tqdm import tqdm                                         # suppress the status bar
from functools import partialmethod
tqdm.__init__ = partialmethod(tqdm.__init__, disable=True)

ignore_warnings = True                                        # ignore warnings?

import matplotlib.pyplot as plt                               # for plotting
import matplotlib as mpl                                      # custom colorbar
import matplotlib.ticker as mticker                           # custom colorpar ticks
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
plt.rc('axes', axisbelow=True)                                # plot all grids below the plot elements
if ignore_warnings == True:
    import warnings
    warnings.filterwarnings('ignore')

import os
import matplotlib.image as mpimg
import subprocess
import time
from matplotlib.offsetbox import AnchoredText
import matplotlib.gridspec as gridspec

import plotly.io as pio
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

# trimmed inverted rainbow colormap hsv removing magenta tones
cmap = mpl.colors.ListedColormap(mpl.cm.hsv(np.linspace(0.0, 0.70, 256))).reversed()

GeostatsPy version: 0.0.78

Nuestro punto de partida sera la base de datos simplificada de la sección anterior. Filtramos aquellas muestras que pertenezcan a nuestra UG de interés. Luego exportamos la base de datos en formato GSLIB para ser procesada por la libreria. La exportación no debe contener ningún caracter y solo debe contener números, sino no podrá ser procesado.

## Cargar data simplificada de la sección anterior

import pandas as pd
import numpy as np

# Cargar base de datos en pandas y aproximar a dos decimales
df = pd.read_csv('data.csv', sep=',', encoding='latin1').round(2)
df['UG'] = df['UG'].astype(int)  # Asegurarse de que la columna 'UG' sea de tipo entero

# Filtrar data para quedarnos solo con UG 4
df = df[df['UG'] == 4].reset_index(drop=True)

# Guardar data con indicadores en formato GSLIB
GSLIB.Dataframe2GSLIB('./01_data/data_UG4.dat',df)  # guardar data

df.head()  # mostrar las primeras filas de la data

	X	Y	Z	cu_pct	UG
0	472204.97	6925835.91	3992.56	0.26	4
1	472205.48	6925836.65	3986.63	0.22	4
2	472205.98	6925837.38	3980.70	0.32	4
3	472206.49	6925838.12	3974.76	0.40	4
4	472206.99	6925838.86	3968.83	0.34	4

# Chequeamos algunas estadísticas básicas de la muestras
df['cu_pct'].describe()

count    8631.000000
mean        0.246391
std         0.222999
min         0.000000
25%         0.120000
50%         0.210000
75%         0.320000
max         6.720000
Name: cu_pct, dtype: float64

4.2. Búsqueda de anisotropías y variograma experimental

4.2.1. Mapas Variográficos

Utilizamos los mapas variográfico para buscar anisotropías. A continuacion se explican los parametros requeridos para correr el programam VARMAP de la GSLIB.

• file with data: archivo con los datos de sondajes.

• number of variables: column numbers: primer numero es el nº de variables a analizar (será 1 pues solo tenemos la ley). Segundo nº indica la columna en el archivo data.dat donde se ubica la variable.

• trimming limits: se indica el rango de valores de la variable que usaremos para hacer el análisis (se usará de preferencia de 0 a infinito).

• 1=regular grid, 0=scattered values: Nos pregunta si los sondajes están en grilla regular ( 1) o dispersos (0).

• if =1: nx, ny, nz: Esta línea se edita solo cuando en la línea anterior se coloca 1, en nuestro caso se deja tal cual. *** xsiz, ysiz, zsiz**: Esta línea se edita solo cuando en la línea anterior se coloca 1, en nuestro caso se deja tal cual.

• if =0: columns for x,y, z coordinates:Se ingresa la columna para la coordenada x,y,z en el archivo de entrada ( Puntosporfido.prn ).

• file for variogram output: Se pone el nombre del archivo de salida de este ejecutable (debe tener extensión *.out).

• nxlag, nylag, nzlag: Se pone el nº de pasos que queremos que tengan los variogramas del mapa variográfico. Se debe recordar que el nº de pasos (nxlag) multiplicado por el largo del paso en una coordenada (dxlag) debe ser igual (aprox) a la mitad de la extensión del dominio en esa coordenada. Por ejemplo, si nuestro deposito mide 2000 m. en coord. X, luego se harán 20 pasos de 50m lo cual es aprox la mitad del largo del dominio.

• dxlag, dylag, dzlag: Serán los largos de los pasos del mapa en las diferentes direcciones.

• minimum number of pairs: Se coloca el mínimo nº de pares de datos que queremos que el mapa requiera para realizar el cálculo en una dirección y distancia (de preferencia se pone 1).

• standardize sill? (0=no, 1=yes): Nos pregunta si queremos una meseta. De preferencia se pone 0.

• number of variograms: Nos pide el número de variogramas a hacer, será 1.

• tail, head, variogram type: Nos pide la cabeza y cola de los datos de variograma en los primeros 2 números, en el 3ro nos pide el tipo de variograma que aparece en la lista de abajo (se debe colocar 1 en todos los ítems)

nxlag=25; nylag=25; nzlag=25
dxlag=50; dylag=50; dzlag=50


parfile = '''                  Parameters for VARMAP
                  *********************

START OF PARAMETERS:
./01_data/data_UG4.dat          -file with data
1   4                        -   number of variables: column numbers
0     1           -   trimming limits
0                            -1=regular grid, 0=scattered values
 50   50    1                -if =1: nx,     ny,   nz - No mover
1.0  1.0  1.0                -       xsiz, ysiz, zsiz  - No mover
1   2   3                    -if =0: columns for x,y, z coordinates
./03_resultados/varmap_UG4.out                   -file for variogram output
'''+str(nxlag)+" "+str(nylag)+" "+str(nzlag)+'''              -nxlag, nylag, nzlag
'''+str(dxlag)+" "+str(dylag)+" "+str(dzlag)+'''              -dxlag, dylag, dzlag
5                            -minimum number of pairs
0                            -standardize sill? (0=no, 1=yes)
1                            -number of variograms
1   1   1                    -tail, head, variogram type
'''

with open("02_parametros/varmap_UG4.par", "w") as f:
    f.write(parfile)

# Ejecutar vmodel de GSLIB
result = subprocess.run(["./Gslib90/varmap", "./02_parametros/varmap_UG4.par"], capture_output=True, text=True, check=True)

El resultado del mapa variográfico tendrá el doble de celdas indicadas en cada dirección 2*nxlag+1, 2*nylag+1, 2*nzlag+1

Usamos los parametros definidos anteriormente para leer la salida del programa VARMAP

nxlag, nylag, nzlag

dxlag, dylag, dzlag

varmap_UG4 = GSLIB.GSLIB2ndarray_3D("./03_resultados/varmap_UG4.out ", 0, 1, 2*nxlag+1, 2*nylag+1, 2*nzlag+1)[0][0]

idz = 25; idy = 25; idx = 25   # índices para cortes en z, y, x respectivamente. Valores corresponden a nxlag/2, nylag/2, nzlag/2

plt.subplots()
plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura | limites escala de colores  ↓   ↓   de 0.0 a 0.03  - Modificar si es necesario
GSLIB.pixelplt_st(varmap_UG4[idz,:,:],-nxlag*dxlag,nxlag*dxlag,-nylag*dylag,nylag*dylag,dzlag,0.0,0.03,'Mapa Variográfico Cu - UG 4 - Vista en planta','X(m)','Y(Mm)','$\mathbf{\gamma}$',cmap)

plt.subplots()
plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura | limites escala de colores  ↓   ↓   de 0.0 a 0.03  - Modificar si es necesario
GSLIB.pixelplt_st(varmap_UG4[:,idy,:],-nylag*dylag,nylag*dylag,-nzlag*dzlag,nzlag*dzlag,dxlag,0.0,0.03,'Mapa Variográfico Cu - UG 4 - Vista Este','Y(m)','Z(m)','$\mathbf{\gamma}$',cmap)

plt.subplots()
plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura | limites escala de colores  ↓   ↓   de 0.0 a 0.03  - Modificar si es necesario
GSLIB.pixelplt_st(varmap_UG4[:,:,idx],-nxlag*dxlag,nxlag*dxlag,-nzlag*dzlag,nzlag*dzlag,dxlag,0.0,0.03,'Mapa Variográfico Cu - UG 4 - Vista Norte','X(m)','Z(m)','$\mathbf{\gamma}$',cmap)

<matplotlib.image.AxesImage at 0x125d30d3110>

4.2.2. Cálculo experimental del variograma

Definimos la siguiente función, que nos permitirá optimizar la ejecución del calculo de variogramas.

def gamv_3d(df, xcol, ycol, zcol, vcol, nlag, tmin,tmax, lagdist,lag_tol, azi, atol, bandh, dip, dtol, bandv, isill):
    """Calcula el variograma experimental 3D usando GSLIB GAMV."""
    lag = []
    gamma = []
    npair = []

    df_ext = pd.DataFrame({"X": df[xcol], "Y": df[ycol], "Z": df[zcol],"Variable": df[vcol]})
    GSLIB.Dataframe2GSLIB("01_data/gamv_out.dat", df_ext)

    with open("02_parametros/gamv.par", "w") as f:
        f.write("                  Parameters for GAMV                                      \n")
        f.write("                  *******************                                      \n")
        f.write("                                                                           \n")
        f.write("START OF PARAMETERS:                                                       \n")
        f.write("01_data/gamv_out.dat                    -file with data                            \n")
        f.write("1   2   3                         -   columns for X, Y, Z coordinates      \n")
        f.write("1   4   0                         -   number of variables,col numbers      \n")
        f.write( str(tmin) + " " + str(tmax)+ "  -   trimming limits                      \n")
        f.write("03_resultados/gamv.out                          -file for variogram output               \n")
        f.write(str(nlag) + "                      -number of lags                          \n")
        f.write(str(lagdist) + "                       -lag separation distance                 \n")
        f.write(str(lag_tol) + "                   -lag tolerance                           \n")
        f.write("1                                 -number of directions                    \n")
        f.write(str(azi) + " " + str(atol) + " " + str(bandh) + " " +str(dip) + " " + str(dtol) + " " + str(bandv)+ "  -azm,atol,bandh,dip,dtol,bandv \n")
        f.write(str(isill) + "                    -standardize sills? (0=no, 1=yes)        \n")
        f.write("1                                 -number of variograms                    \n")
        f.write("1   1   1                         -tail var., head var., variogram type    \n")

    # Ejecutar gamv de GSLIB
    result = subprocess.run(["./Gslib90/gamv", "./02_parametros/gamv.par"], capture_output=True, text=True, check=True)

    with open("03_resultados/gamv.out") as f:
        next(f)  # skip the first line

        for line in f:
            _, l, g, n, *_ = line.split()
            lag.append(float(l))
            gamma.append(float(g))
            npair.append(float(n))

    return lag, gamma, npair, result

def make_variogram3D(
    nug,
    nst,
    it1,
    cc1,
    azi1,
    dip1,
    hmaj1,
    hmin1,
    hvert1,
    it2=1,
    cc2=0,
    azi2=0,
    dip2=0,
    hmaj2=0,
    hmin2=0,
    hvert2=0,
    it3=1,
    cc3=0,
    azi3=0,
    dip3=0,
    hmaj3=0,
    hmin3=0,
    hvert3=0
):
    """Crea un variograma 3D con tres estructuras.
    """

    if cc2 == 0:
        nst = 1
    var = dict(
        [
            ("nug", nug),
            ("nst", nst),
            ("it", np.array([it1,it2,it3])),
            ("cc", np.array([cc1,cc2,cc3])),
            ("azi", np.array([azi1,azi2,azi3])),
            ("dip", np.array([dip1,dip2,dip3])),
            ("hmaj", np.array([hmaj1,hmaj2,hmaj3])),
            ("hmin", np.array([hmin1,hmin2,hmin3])),
            ("hvert", np.array([hvert1,hvert2,hvert3])),
        ]
    )

    return var

def vmodel(nlag,xlag,mazm,mdip,vario):
    """Calcula el variograma experimental 3D usando GSLIB VMODEL."""
    lag = []
    gamma = []

    with open("02_parametros/vmodel.par", "w") as f:
        f.write("                  Parameters for VMODEL                                      \n")
        f.write("                  *******************                                      \n")
        f.write("                                                                           \n")
        f.write("START OF PARAMETERS:                                                     \n")
        f.write("./03_resultados/vmodel.out                   -file for variogram output                 \n")
        f.write("1 " + str(nlag)+ "                          -number of directions and lags             \n")
        f.write(str(mazm) + " " + str(mdip) + " " + str(xlag) + "  -azm, dip, lag distance                     \n")
        f.write(str(vario["nst"]) + " " + str(vario['nug']) + "  -nst, nugget effect                        \n")
        f.write(str(vario["it"][0]) + " " + str(vario["cc"][0]) + " " + str(vario["azi"][0]) + " " + str(vario["dip"][0])  + "  0.0    -it,cc,ang1,ang2,ang3                       \n")
        f.write("         " + str(vario["hmaj"][0]) + " " + str(vario["hmin"][0]) + " " + str(vario["hvert"][0])  + "  -a_hmax, a_hmin, a_vert                    \n")
        f.write(str(vario["it"][1]) + " " + str(vario["cc"][1]) + " " + str(vario["azi"][1]) + " " + str(vario["dip"][1]) + "  0.0    -it,cc,ang1,ang2,ang3                       \n")
        f.write("         " + str(vario["hmaj"][1]) + " " + str(vario["hmin"][1]) + " " + str(vario["hvert"][1])  + "  -a_hmax, a_hmin, a_vert                    \n")
        f.write(str(vario["it"][2]) + " " + str(vario["cc"][2]) + " " + str(vario["azi"][2]) + " " + str(vario["dip"][2]) + "  0.0    -it,cc,ang1,ang2,ang3                       \n")
        f.write("         " + str(vario["hmaj"][2]) + " " + str(vario["hmin"][2]) + " " + str(vario["hvert"][2])  + "  -a_hmax, a_hmin, a_vert                    \n")

    # Ejecutar vmodel de GSLIB
    result = subprocess.run(["./Gslib90/vmodel", "./02_parametros/vmodel.par"], capture_output=True, text=True, check=True)

    with open("03_resultados/vmodel.out") as f:
        next(f)  # skip the first line

        for line in f:
            _, l, g, _, _, *_ = line.split()
            lag.append(float(l))
            gamma.append(float(g))

    return lag, gamma, result

Estamos listos para calcular los variogramas. Vamos a calcular variogramas isótropos para todas las UGs. Algunas notas sobre los parámetros:

• number of variables,col numbers indica el numero de variables a considerar (1, solo la variable de cobre en nuestro caso) y en que columna se encuentra esta variable en el archivo

• trimming limits son límites de recorte — establecidos para no tener impacto, no es necesario filtrar los datos.

• lag separation distance, lag tolerance son la distancia de lag y la tolerancia de lag — definidos según el espaciado común de los datos (50 m) y con una tolerancia del 100% de la distancia de lag para un mayor suavizado.

• number of lags es el número de lags — definido para extenderse un poco más allá del 20% del rango de los datos.

• bandh es el ancho de banda horizontal — configurado para no tener efecto.

• bandv es el ancho de banda vertical — configurado para no tener efecto.

• azi es el acimut — no tiene efecto porque hemos fijado atol (tolerancia de acimut) en 90.0.

• standardize sills? es un indicador para estandarizar la distribución a varianza 1 — no será usado durante este curso.

• number of variograms indica el numero de variogramas que calcularemos. 3 en nuestro caso (uno para cada indicador).

• tail var., head var., variogram type finalmente indicamos que variables entran en el cálculo del variograma y el tipo (variograma tradicional si el indicador es 1)

Probemos los siguiente parámetros para el cálculo del variograma en distintas direcciones (buscando alguna anisotropía) y visualicemos los resultados.

tmin = 0.; tmax = 1.5;
lag_dist = 30.0; lag_tol = 15.0; nlag = 10; bandh = 60; atol = 22.5; isill = 0;
dip = 60; dtol = 12.5; bandv = 20; ## Mover estos parámetros a gusto del usuario

A continuación, solo modificamos la variable azi

# Variograma experimental y búsqueda de anisotropías entre 0 y 165 grados

continuous = "cu_pct" # nombre de la variable continua

lag_0, gamma_0, npair_0, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 0, atol, bandh, dip, dtol, bandv, isill)
lag_15, gamma_15, npair_15, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 15, atol, bandh, dip, dtol, bandv, isill)
lag_30, gamma_30, npair_30, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 30, atol, bandh, dip, dtol, bandv, isill)
lag_45, gamma_45, npair_45, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 45, atol, bandh, dip, dtol, bandv, isill)
lag_60, gamma_60, npair_60, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 60, atol, bandh, dip, dtol, bandv, isill)
lag_75, gamma_75, npair_75, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 75, atol, bandh, dip, dtol, bandv, isill)
lag_90, gamma_90, npair_90, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 90, atol, bandh, dip, dtol, bandv, isill)
lag_105, gamma_105, npair_105, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 105, atol, bandh, dip, dtol, bandv, isill)
lag_120, gamma_120, npair_120, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 120, atol, bandh, dip, dtol, bandv, isill)
lag_135, gamma_135, npair_135, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 135, atol, bandh, dip, dtol, bandv, isill)
lag_150, gamma_150, npair_150, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 150, atol, bandh, dip, dtol, bandv, isill)
lag_165, gamma_165, npair_165, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 165, atol, bandh, dip, dtol, bandv, isill)

fig, ax = plt.subplots(figsize=(5,2.5))
ax.clear()

plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
plt.rcParams["font.family"] = "Times New Roman"

plt.plot(lag_0,gamma_0,color = 'red',label = 'azi = 0°',zorder=8)
plt.plot(lag_15,gamma_15,color = 'orange',label = 'azi = 15°',zorder=8)
plt.plot(lag_30,gamma_30,color = 'yellow',label = 'azi = 30°',zorder=8)
plt.plot(lag_45,gamma_45,color = 'green',label = 'azi = 45°',zorder=8)
plt.plot(lag_60,gamma_60,color = 'blue',label = 'azi = 60°',zorder=8)
plt.plot(lag_75,gamma_75,color = 'indigo',label = 'azi = 75°',zorder=8)
plt.plot(lag_90,gamma_90,color = 'violet',label = 'azi = 90°',zorder=8)
plt.plot(lag_105,gamma_105,color = 'brown',label = 'azi = 105°',zorder=8)
plt.plot(lag_120,gamma_120,color = 'pink',label = 'azi = 120°',zorder=8)
plt.plot(lag_135,gamma_135,color = 'gray',label = 'azi = 135°',zorder=8)
plt.plot(lag_150,gamma_150,color = 'cyan',label = 'azi = 150°',zorder=8)
plt.plot(lag_165,gamma_165,color = 'magenta',label = 'azi = 165°',zorder=8)

plt.xlabel(r'Distancia (m), $\bf(h)$ '); plt.ylabel(r'$\gamma \bf(h)$')
plt.title('Variograma de ley de Cu en UG 4')
plt.xlim([0,lag_dist*nlag]); plt.ylim([0,0.04]); plt.legend(loc='upper left');

## add grid lines
plt.grid(visible=True, which='major', color='lightgrey', linestyle='-', linewidth=0.5)
plt.grid(visible=True, which='minor', color='lightgrey', linestyle='--', linewidth=0.2)
plt.minorticks_on()

# mover legenda fuera del plot
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)


plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.3)
plt.show()

Estos variogramas son bastante feos pero igualmente nos permiten identificar la dirección de mayor continuidad, la cual identificamos como la dirección N45E, con mayor alcance.

Reajustamos nuestros parámetros y calculamos esta vez el variograma en esta orientación y su orientación perpendicular N135E.

lag_dist = 30.0; lag_tol = 15; nlag = 10; bandh = 50; atol = 22.5; isill = 0;
dip = 60; dtol = 12.5; bandv = 50; ## Mover estos parámetros a gusto del usuario

continuous = "cu_pct" # nombre de la variable continua

az_max = 45 # dirección de máxima continuidad
az_min = 135  # dirección de mínima continuidad

lag_max, gamma_max, npair_max, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, az_max, atol, bandh, dip, dtol, bandv, isill)
lag_min, gamma_min, npair_min, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, az_min, atol, bandh, dip, dtol, bandv, isill)

fig, ax = plt.subplots(figsize=(5,2.5))
ax.clear()

plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
plt.rcParams["font.family"] = "Times New Roman"


plt.plot(lag_max[1:],gamma_max[1:],color = 'red',label = 'azi = '+str(az_max)+'°',zorder=9)
plt.plot(lag_min[1:],gamma_min[1:],color = 'blue',label = 'azi = '+str(az_min)+'°',zorder=9)

plt.scatter(lag_max[1:],gamma_max[1:],color = 'red',edgecolor='black',s=np.multiply(npair_max[1:],0.01),marker='o',zorder=8)
plt.scatter(lag_min[1:],gamma_min[1:],color = 'blue',edgecolor='black',s=np.multiply(npair_min[1:],0.01),marker='o',zorder=8)

plt.xlabel(r'Distancia (m), $\bf(h)$ '); plt.ylabel(r'$\gamma \bf(h)$')
plt.title('Variograma de ley de Cu en UG 4')
plt.xlim([0,lag_dist*nlag]); plt.ylim([0,0.04]); plt.legend(loc='upper left');

## add grid lines
plt.grid(visible=True, which='major', color='lightgrey', linestyle='-', linewidth=0.5)
plt.grid(visible=True, which='minor', color='lightgrey', linestyle='--', linewidth=0.2)
plt.minorticks_on()

# mover legenda fuera del plot
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)


plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.3)
plt.show()

Continuamos calculando el variograma en la dirección vertical.

lag_dist = 20.0; lag_tol = 20.0; nlag = 15; bandh = 30; atol = 180; isill = 0;
dip_vert = 150; dtol = 15; bandv = 30; ## Mover estos parámetros a gusto del usuario

# Variograma experimental y búsqueda de anisotropías entre 0 y 165 grados

continuous = "cu_pct" # nombre de la variable continua

lag_vert, gamma_vert, npair_vert, result = gamv_3d(df,"X","Y","Z",continuous, nlag, tmin,tmax, lag_dist,lag_tol, 0, atol, bandh, dip_vert, dtol, bandv, isill)

fig, ax = plt.subplots(figsize=(5,2.5))
ax.clear()

plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
plt.rcParams["font.family"] = "Times New Roman"

plt.plot(lag_max[1:],gamma_max[1:],color = 'red',label = 'azi = '+str(az_max)+'°',zorder=9)
plt.plot(lag_min[1:],gamma_min[1:],color = 'blue',label = 'azi = '+str(az_min)+'°',zorder=9)
plt.plot(lag_vert[1:],gamma_vert[1:],color = 'black',label = 'dip = '+str(dip_vert)+'°',zorder=9)


plt.scatter(lag_max[1:],gamma_max[1:],color = 'red',edgecolor='black',s=np.multiply(npair_max[1:],0.01),marker='o',zorder=8)
plt.scatter(lag_min[1:],gamma_min[1:],color = 'blue',edgecolor='black',s=np.multiply(npair_min[1:],0.01),marker='o',zorder=8)
plt.scatter(lag_vert[1:],gamma_vert[1:],color = 'black',edgecolor='black',s=np.multiply(npair_vert[1:],0.01),marker='o',zorder=8)

plt.xlabel(r'Distancia (m), $\bf{h}$ '); plt.ylabel(r'$\gamma \bf(h)$')
plt.title('Variograma Cu - UG 4')
plt.xlim([0,lag_dist*nlag]); plt.ylim([0,0.04]); plt.legend(loc='upper left');

## add grid lines
plt.grid(visible=True, which='major', color='lightgrey', linestyle='-', linewidth=0.5)
plt.grid(visible=True, which='minor', color='lightgrey', linestyle='--', linewidth=0.2)
plt.minorticks_on()

# mover legenda fuera del plot
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)


plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.3)
plt.show()

Usamos la función make_variogram de GeostatsPy para construir un objeto modelo de variograma.

• Es un diccionario para almacenar de forma compacta los parámetros del modelo de variograma que se usan en la visualización, kriging y simulación.

Los parámetros del modelo de variograma incluyen:

• nug - aporte del efecto pepita (nugget) a la meseta (sill)

• nst - número de estructuras anidadas (1 o 2)

• it - tipo para la estructura (1 - esférico, 2 - exponencial, 3 - gaussiano)

• cc - contribución de cada estructura anidada (las contribuciones + nugget deben sumar la meseta)

• azi - acimut para esta estructura (dirección mayor), la minor es ortogonal

• hmaj - alcance en la dirección mayor

• hmin - alcance en la dirección menor

Usamos arreglos para it, cc, azi, hmaj y hmin cuando hay más de una estructura anidada.

• Si solo usamos 1 estructura (más opción de nugget), omita los parámetros de la segunda estructura y estos por defecto serán \(cc2 = 0\) (sin contribución).

Tomemos como ejemplo la UG 2. Aquí está mi modelo:

nug = 0.003; nst = 3     # definimos el efecto pepa (0.0) y el número de estructuras anidadas (1)
it1 = 2; cc1 = 0.022; azi1 = az_max; dip1=dip; hmaj1 = 200; hmin1 = 100; hvert1 = 90;  # 1 estructura anidada (exponencial)
it2 = 2; cc2 = 0.005; azi2 = az_max; dip2=dip; hmaj2 = 9999; hmin2 = 100; hvert2 = 90;  # 1 estructura anidada (exponencial)
it3 = 2; cc3 = 0.005; azi3 = az_max; dip3=dip; hmaj3 = 9999; hmin3 = 9999; hvert3 = 90;  # 1 estructura anidada (exponencial)

vario_cu_UG4 = make_variogram3D(nug,nst,it1,cc1,azi1,dip1,hmaj1,hmin1,hvert1,it2,cc2,azi2,dip2,hmaj2,hmin2,hvert2,it3,cc3,azi3,dip3,hmaj3,hmin3,hvert3)  # crear el objeto del modelo de variograma

4.2.3. Trazado del modelo de variograma

Para trazar el variograma usamos el programa vmodel que proyecta el modelo a un conjunto de distancias de lag en las direcciones mayor y menor.

Las entradas para vmodel son:

• nlag - número de puntos a lo largo del variograma para la proyección

• xlag - tamaño del lag para la proyección

• azm - dirección de la proyección en azimut (trabajamos en 2D)

• vario - el diccionario del modelo de variograma creado con make_variogram

Nota: esta función es solo para visualización; por eso es habitual usar un xlag muy pequeño y un nlag grande para una visualización de alta resolución del modelo.

Las salidas de vmodel incluyen:

• lag distance - distancia offset a lo largo de la proyección (el \(\mathbf{h}\) en el gráfico)

• variogram - valor del variograma en esa distancia de lag (el \(\gamma(\mathbf{h})\))

• resutls - verificación de la correcta ejecución (imprimir en caso de que algo ande mal)

Tenemos 2 estructuras y sin efecto pepa como ejemplo para el caso de la UG 2.

nlag = 50; xlag = 5; a_max = az_max; dip_p = dip;  # project the model in max range azimuth
h_max , gam_max , result = vmodel(nlag,xlag,a_max,dip_p,vario_cu_UG4)

nlag = 50; xlag = 5; a_min = az_min; dip_p = 0;  # project the model in min range azimuth
h_min , gam_min , result = vmodel(nlag,xlag,a_min,dip_p,vario_cu_UG4)

nlag = 50; xlag = 5; a_vert = az_max; dip_p = dip_vert;  # project the model in vertical direction
h_vert , gam_vert , result = vmodel(nlag,xlag,a_vert,dip_p,vario_cu_UG4)

fig, ax = plt.subplots(figsize=(5,2.5))
ax.clear()

plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
plt.rcParams["font.family"] = "Times New Roman"

plt.scatter(lag_max[1:],gamma_max[1:],color = 'red',edgecolor='black',s=np.multiply(npair_max[1:],0.02),marker='o',zorder=8)
plt.scatter(lag_min[1:],gamma_min[1:],color = 'blue',edgecolor='black',s=np.multiply(npair_min[1:],0.02),marker='o',zorder=8)
plt.scatter(lag_vert[1:],gamma_vert[1:],color = 'black',edgecolor='black',s=np.multiply(npair_vert[1:],0.01),marker='o',zorder=8)

plt.plot(h_max,gam_max,color = 'red',label = 'ajuste azi = '+str(az_max)+'°',zorder=9)
plt.plot(h_min,gam_min,color = 'blue',label = 'ajuste azi = '+str(az_min)+'°',zorder=9)
plt.plot(h_vert,gam_vert,color = 'black',label = 'ajuste vertical, dip = '+str(dip_p)+'°',zorder=9)

plt.xlabel(r'Distancia (m), $\bf(h)$ '); plt.ylabel(r'$\gamma \bf(h)$')
plt.title('Modelo Variográfico de Cu - UG 4')
plt.xlim([0,nlag*xlag]); plt.ylim([0,0.045]); plt.legend(loc='upper left');

## add grid lines
plt.grid(visible=True, which='major', color='lightgrey', linestyle='-', linewidth=0.5)
plt.grid(visible=True, which='minor', color='lightgrey', linestyle='--', linewidth=0.2)
plt.minorticks_on()

# mover legenda fuera del plot
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.2, wspace=0.2, hspace=0.3)
plt.show()

4.3. Selección de plan de estimación

4.3.1. Validación cruzada mediante Kriging para variables continuas - Ley de Cu - UG 4

Aplicamos los conceptos desarrollados anteriormente y el modelo de variograma obtenido para cada la ley de Cu.

Este capítulo presenta los programas de kriging de GSLIB. Estos programas permiten la validación cruzada mediante kriging simple (SK), kriging ordinario (OK) y kriging con varios modelos de tendencia (KT).

El programa kt3d ofrece un programa de kriging 3D bastante avanzado para puntos o bloques mediante kriging simple (SK), kriging ordinario (OK) o kriging con un modelo de tendencia polinómica (no sera revisado en este curso).

Los parámetros requeridos por el programa principal kt3d se listan a continuación:

• datafl: archivo de entrada con los datos en formato GSLIB.

• DH, X, Y, Z, icolvr, icolsec: columnas para la columna que contiene los codigos DH (ya los borramos y no serán necesarios - valor 0) coordenadas x, y, z, la variable a estimar y la variable de deriva externa (o media no estacionaria - no utilizada en este curso).

• tmin y tmax: triming limits o valores fuera de este rango (menores estrictamente que tmin o mayores o iguales que tmax) se ignoran.: valores fuera de este rango (menores estrictamente que tmin o mayores o iguales que tmax) se ignoran..

• option: 0 para kriging sobre una rejilla de puntos o bloques, 1 para valida validación cruzada con los datos en datafl, 2 para jackknife con los datos del archivo siguiente.

• jackfl: archivo con las ubicaciones para realizar estimaciones (opción jackknife).

• icoljx, icoljy, icoljz, icoljvr, icoljsec: columnas para x, y, z, la variable y la variable secundaria en jackfl.kfl.

• idbg: nivel entero de depuración entre 0 y 3. A mayor nivel, más información entregada. Los niveles usuales son 0 y 1; 2 y 3 imprimen todas las matrices de kriging y datos usados para cada punto/bloque (no recomendable con grillas grandes).

• dbgfl: archivo donde se escribe la salida de depuración.

• outfl: archivo de salida de la rejilla. Contiene la estimación y la varianza de kriging para cada punto/bloque en la rejilla, con orden rápido en x, luego y y finalmente z. Los puntos no estimados se marcan con un número negativo grande (-999.).

• nx, xmn, xsiz: definición de la grilla en la dirección x.

• ny, ymn, ysiz: definición de la grilla en la dirección y.

• nz, zmn, zsiz: definición de la grilla en la dirección z.

• nxdis, nydis, nzdis: número de puntos de discretización para un bloque. Si nxdis, nydis y nzdis son 1, se realiza kriging puntual.

• ndmin y ndmax: número mínimo y máximo de datos a usar para kriging de un bloque.

• noct: máximo a retener por octante (no se usa búsqueda por octantes si noct = 0).i noct = 0).

• radius_hmax, radius_hmin, radius_vert: radios de búsqueda en la dirección horizontal máxima, horizontal mínima y vertical.

• sang1, sang2, sang3: ángulos que describen la orientación de la elipsoide de búsqueda (anisotropía).

• ikrige y skmean: si ikrige = 0 se realiza kriging simple estacionario con media skmean; si ikrige = 1 se realiza kriging ordinario y el valor de skmean no es considerado.

• nst y c0: número de estructuras del variograma y constante nugget. El efecto pepita no se cuenta como estructura.

Para cada una de las nst estructuras anidadas se debe definir el tipo (it), la contribución (cc), los ángulos (azi, dip) y los alcances (a_hmax, a_hmin, a_vert).

4.3.1.1. Plan de búsqueda sin uso de octantes

Vemos como queda nuestro archivo de parámetros para la UG 4, en un primer caso sin consideración de octantes, y buscando 24 datos.

parfile = '''                  Parameters for KT3D
                  *******************

START OF PARAMETERS:
./01_data/data_UG4.dat              -file with data
0  1  2  3  4  0                 -   columns for DH,X,Y,Z,var,sec var
-1.0e21   1.0e21                 -   trimming limits
1                                -option: 0=grid, 1=cross, 2=jackknife - AR : usamos la opción 1 de cross-validation
xvk.dat                          -file with jackknife data - AR : no se usa en este caso
1   2   0    3    0              -   columns for X,Y,Z,vr and sec var - AR : no se usa en este caso
0                                -debugging level: 0,1,2,3
./03_resultados/kt3d-pt-cross_UG4.dbg                   -file for debugging output
./03_resultados/kt3d-pt-cross_UG4.out                   -file for kriged output
128   471116.24    15                       -nx,xmn,xsiz - AR : no se usa en este caso debido a cross-validation
97   6924721.57    15                       -ny,ymn,ysiz - AR : no se usa en este caso debido a cross-validation
100   2950.02    15                     -nz,zmn,zsiz - AR : no se usa en este caso debido a cross-validation
1    1      1                    -x,y and z block discretization
20    24                           -min, max data for kriging - AR : uso de 8 datos para el krigeo
0                                -max per octant (0-> not used) - AR : 0 si no se desea usar octantes
200.0  100.0  90.0                -maximum search radii - AR : radios de búsqueda de acuerdo a alcances de variograma
 45   60   0.0                 -angles for search ellipsoid - AR : ángulos de búsqueda de acuerdo a azimuts y dips principales del variograma
1     0.37                      -0=SK,1=OK,2=non-st SK,3=exdrift - AR : uso de kriging ordinario, segundo valor no es relevante en este caso
0 0 0 0 0 0 0 0 0                -drift: x,y,z,xx,yy,zz,xy,xz,zy
0                                -0, variable; 1, estimate trend
extdrift.dat                     -gridded file with drift/mean
4                                -  column number in gridded file
3 0.003  -nst, nugget effect  - AR: variograma de CU para UG 4  modelado en sección anterior
2 0.022 45 60  0.0    -it,cc,ang1,ang2,ang3
         200 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 9999 90  -a_hmax, a_hmin, a_vert
'''

with open("02_parametros/kt3d_cross_UG4.par", "w") as f:
    f.write(parfile)

Ahora ya podemos ejecutar kt3d de GSLIB. Puede tardar varios minutos dependiendo del tamaño de la grilla definida.

# Ejecutar kt3d de GSLIB
result = subprocess.run(["./Gslib90/kt3d", "./02_parametros/kt3d_cross_UG4.par"], capture_output=True, text=True, check=True)

ValReal = []
Estim = []
Error = []

with open("./03_resultados/kt3d-pt-cross_UG4.out") as f:
    head = [next(f) for _ in range(2)]  # read first two lines
    line2 = head[1].split()
    ncol = int(line2[0])  # get the number of columns

    for icol in range(ncol):  # read over the column names
        head = next(f)

    for line in f:
        _,_,_, l, m, _, n = line.split()
        # skip si es -999
        if float(m) < -99:
            continue
        ValReal.append(float(l))
        Estim.append(float(m))
        Error.append(float(n))

## Plot de validación cruzada
import plotly.graph_objects as go

# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode

fig = go.Figure()
fig.add_trace(go.Scatter(x=ValReal, y=Estim, mode='markers',
                         marker=dict(color='black', size=5, opacity=0.7),
                         name='Puntos de Validación Cruzada'))

# Añadir línea de 1:1 entre 0 y maximo valor de ValReal y Estim
max_val = max(max(ValReal), max(Estim))
fig.add_trace(go.Scatter(x=[0, max_val], y=[0, max_val], mode='lines',
                         line=dict(color='red', dash='dash'),
                         name='Línea 1:1'))

# Agregar información estadística en recuadro hacia abajo
media_real = sum(ValReal) / len(ValReal)
media_estim = sum(Estim) / len(Estim)
var_real = sum((x - media_real) ** 2 for x in ValReal) / len(ValReal)
var_estim = sum((x - media_estim) ** 2 for x in Estim) / len(Estim)
corr = np.corrcoef(ValReal, Estim)[0, 1]

stats_text = (f'Numero de Muestras: {len(ValReal)}\n'
              f'Media Real: {media_real:.2f}\n'
              f'Media Estimada: {media_estim:.2f}\n'
              f'Varianza Real: {var_real:.2f}\n'
              f'Varianza Estimada: {var_estim:.2f}\n'
              f'Correlación: {corr:.2f}\n'
              f'Rango Real: {min(ValReal):.2f} - {max(ValReal):.2f}\n'
              f'Rango Estimado: {min(Estim):.2f} - {max(Estim):.2f}'
              )

fig.add_annotation(
    x=max_val * 0.85,
    y=max_val * 0.25,
    xref="x",
    yref="y",
    text=stats_text.replace('\n', '<br>'),
    showarrow=False,
    align="left",
    font=dict(family="Times New Roman", size=20, color="black"),
    bordercolor="black",
    borderwidth=1,
    borderpad=10,
    bgcolor="white",
    opacity=0.8
)

fig.update_layout(
    title='Val. Cruzada K.O. - Cu - UG 4 - 24 muestras, sin Oct.',
    xaxis_title='Valor Real (Cu %)',
    yaxis_title='Valor Estimado (Cu %)',
    xaxis=dict(range=[0, max_val], scaleanchor="y", scaleratio=1),
    yaxis=dict(range=[0, max_val]),
    width=600,
    height=600,
    showlegend=False,
    hovermode=False,
    font=dict(family="Times New Roman", size=18, color="black"),
    plot_bgcolor='white',
    xaxis_gridcolor='lightgrey',
    yaxis_gridcolor='lightgrey'
)
fig.show()

# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
# Histograma de errores
media_error = sum(Error) / len(Error)
var_error = sum((x - media_error) ** 2 for x in Error) / len(Error)

# Histograma de errores (eje x entre percentil 1 y 99), 10 barras, con recuadro de estadísticas
p1, p99 = np.percentile(Error, [1, 99])
errs = np.array(Error)
errs_clip = errs[(errs >= p1) & (errs <= p99)]

fig = go.Figure()
fig.add_trace(go.Histogram(
    x=errs_clip,
    nbinsx=10,
    marker=dict(color='black', line=dict(color='black', width=1))
))

stats_text = f"Media Error: {media_error:.2f}<br>Var Error: {var_error:.2f}<br>Rango: {min(Error):.2f} , {max(Error):.2f}"

fig.add_annotation(
    x=1.1, y=0.95, xref="paper", yref="paper",
    text=stats_text, showarrow=False, align="left",
    font=dict(family="Times New Roman", size=15, color="black"),
    bordercolor="black", borderwidth=1, borderpad=8, bgcolor="white", opacity=0.9
)

fig.update_layout(
    title='His. de Error - K.O. - Cu - UG 4 - 24 muestras, sin Oct.',
    xaxis_title='Error de Estimación (Cu %)',
    yaxis_title='Frecuencia',
    xaxis=dict(range=[p1, p99], gridcolor='lightgrey'),
    yaxis=dict(gridcolor='lightgrey'),
    plot_bgcolor='white',
    width=500, height=400,
    font=dict(family="Times New Roman", size=15, color="black"),
    showlegend=False
)

fig.show()

4.3.1.2. Plan de búsqueda con uso de octantes

Vemos como queda nuestro archivo de parámetros para la UG 4, en un segundo caso considerando el uso de octantes, y buscando 24 datos.

parfile = '''                  Parameters for KT3D
                  *******************

START OF PARAMETERS:
./01_data/data_UG4.dat              -file with data
0  1  2  3  4  0                 -   columns for DH,X,Y,Z,var,sec var
-1.0e21   1.0e21                 -   trimming limits
1                                -option: 0=grid, 1=cross, 2=jackknife - AR : usamos la opción 1 de cross-validation
xvk.dat                          -file with jackknife data - AR : no se usa en este caso
1   2   0    3    0              -   columns for X,Y,Z,vr and sec var - AR : no se usa en este caso
0                                -debugging level: 0,1,2,3
./03_resultados/kt3d-pt-cross_UG4.dbg                   -file for debugging output
./03_resultados/kt3d-pt-cross_UG4.out                   -file for kriged output
128   471116.24    15                       -nx,xmn,xsiz - AR : no se usa en este caso debido a cross-validation
97   6924721.57    15                       -ny,ymn,ysiz - AR : no se usa en este caso debido a cross-validation
100   2950.02    15                     -nz,zmn,zsiz - AR : no se usa en este caso debido a cross-validation
1    1      1                    -x,y and z block discretization
20    24                           -min, max data for kriging - AR : uso de 24 datos para el krigeo
3                                -max per octant (0-> not used) - AR : número de muestras por octante
200.0  100.0  90.0                -maximum search radii - AR : radios de búsqueda de acuerdo a alcances de variograma
 45   60   0.0                 -angles for search ellipsoid - AR : ángulos de búsqueda de acuerdo a azimuts y dips principales del variograma
1     0.37                      -0=SK,1=OK,2=non-st SK,3=exdrift - AR : uso de kriging ordinario, segundo valor no es relevante en este caso
0 0 0 0 0 0 0 0 0                -drift: x,y,z,xx,yy,zz,xy,xz,zy
0                                -0, variable; 1, estimate trend
extdrift.dat                     -gridded file with drift/mean
4                                -  column number in gridded file
3 0.003  -nst, nugget effect  - AR: variograma de CU para UG 4  modelado en sección anterior
2 0.022 45 60  0.0    -it,cc,ang1,ang2,ang3
         200 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 9999 90  -a_hmax, a_hmin, a_vert
'''

with open("02_parametros/kt3d_cross_UG4.par", "w") as f:
    f.write(parfile)

Ahora ya podemos ejecutar kt3d de GSLIB. Puede tardar varios minutos dependiendo del tamaño de la grilla definida.

# Ejecutar kt3d de GSLIB
result = subprocess.run(["./Gslib90/kt3d", "./02_parametros/kt3d_cross_UG4.par"], capture_output=True, text=True, check=True)

ValReal = []
Estim = []
Error = []

with open("./03_resultados/kt3d-pt-cross_UG4.out") as f:
    head = [next(f) for _ in range(2)]  # read first two lines
    line2 = head[1].split()
    ncol = int(line2[0])  # get the number of columns

    for icol in range(ncol):  # read over the column names
        head = next(f)

    for line in f:
        _,_,_, l, m, _, n = line.split()
        # skip si es -999
        if float(m) < -99:
            continue
        ValReal.append(float(l))
        Estim.append(float(m))
        Error.append(float(n))

## Plot de validación cruzada
import plotly.graph_objects as go

# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode

fig = go.Figure()
fig.add_trace(go.Scatter(x=ValReal, y=Estim, mode='markers',
                         marker=dict(color='black', size=5, opacity=0.7),
                         name='Puntos de Validación Cruzada'))

# Añadir línea de 1:1 entre 0 y maximo valor de ValReal y Estim
max_val = max(max(ValReal), max(Estim))
fig.add_trace(go.Scatter(x=[0, max_val], y=[0, max_val], mode='lines',
                         line=dict(color='red', dash='dash'),
                         name='Línea 1:1'))

# Agregar información estadística en recuadro hacia abajo
media_real = sum(ValReal) / len(ValReal)
media_estim = sum(Estim) / len(Estim)
var_real = sum((x - media_real) ** 2 for x in ValReal) / len(ValReal)
var_estim = sum((x - media_estim) ** 2 for x in Estim) / len(Estim)
corr = np.corrcoef(ValReal, Estim)[0, 1]

stats_text = (f'Numero de Muestras: {len(ValReal)}\n'
              f'Media Real: {media_real:.2f}\n'
              f'Media Estimada: {media_estim:.2f}\n'
              f'Varianza Real: {var_real:.2f}\n'
              f'Varianza Estimada: {var_estim:.2f}\n'
              f'Correlación: {corr:.2f}\n'
              f'Rango Real: {min(ValReal):.2f} - {max(ValReal):.2f}\n'
              f'Rango Estimado: {min(Estim):.2f} - {max(Estim):.2f}'
              )

fig.add_annotation(
    x=max_val * 0.85,
    y=max_val * 0.25,
    xref="x",
    yref="y",
    text=stats_text.replace('\n', '<br>'),
    showarrow=False,
    align="left",
    font=dict(family="Times New Roman", size=20, color="black"),
    bordercolor="black",
    borderwidth=1,
    borderpad=10,
    bgcolor="white",
    opacity=0.8
)

fig.update_layout(
    title='Val. Cruzada K.O. - Cu - UG 4 - 24 muestras, con Oct.',
    xaxis_title='Valor Real (Cu %)',
    yaxis_title='Valor Estimado (Cu %)',
    xaxis=dict(range=[0, max_val], scaleanchor="y", scaleratio=1),
    yaxis=dict(range=[0, max_val]),
    width=600,
    height=600,
    showlegend=False,
    hovermode=False,
    font=dict(family="Times New Roman", size=18, color="black"),
    plot_bgcolor='white',
    xaxis_gridcolor='lightgrey',
    yaxis_gridcolor='lightgrey'
)
fig.show()

# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
# Histograma de errores
media_error = sum(Error) / len(Error)
var_error = sum((x - media_error) ** 2 for x in Error) / len(Error)

# Histograma de errores (eje x entre percentil 1 y 99), 10 barras, con recuadro de estadísticas
p1, p99 = np.percentile(Error, [1, 99])
errs = np.array(Error)
errs_clip = errs[(errs >= p1) & (errs <= p99)]

fig = go.Figure()
fig.add_trace(go.Histogram(
    x=errs_clip,
    nbinsx=10,
    marker=dict(color='black', line=dict(color='black', width=1))
))

stats_text = f"Media Error: {media_error:.2f}<br>Var Error: {var_error:.2f}<br>Rango: {min(Error):.2f} , {max(Error):.2f}"

fig.add_annotation(
    x=1.1, y=0.95, xref="paper", yref="paper",
    text=stats_text, showarrow=False, align="left",
    font=dict(family="Times New Roman", size=15, color="black"),
    bordercolor="black", borderwidth=1, borderpad=8, bgcolor="white", opacity=0.9
)

fig.update_layout(
    title='His. de Error - K.O. - Cu - UG 4 - 24 muestras, con Oct.',
    xaxis_title='Error de Estimación (Cu %)',
    yaxis_title='Frecuencia',
    xaxis=dict(range=[p1, p99], gridcolor='lightgrey'),
    yaxis=dict(gridcolor='lightgrey'),
    plot_bgcolor='white',
    width=500, height=400,
    font=dict(family="Times New Roman", size=15, color="black"),
    showlegend=False
)

fig.show()

4.4. Parámetros de la malla 3D de estimación

Antes de continuar con la estimación, definamos los parámetros de la malla 3D donde se realizará la estimación de las UGs, calculadas ya en la sección anterior.

xmin = 471116.24; xmax = 473023.89                  # rango de valores x
ymin = 6924721.57; ymax = 6926164.63                        # rango de valores y
zmin = 2950.02; zmax = 4445.78                              # rango de valores z

xsiz = 15; ysiz = 15; zsiz = 15                             # tamaño de celda
nx = 128; ny = 97; nz = 100                         # número de celdas

tmin = -999; tmax = 999;                                      # límites de recorte de datos

icx = int(nx/2); cx = xmin + icx*xsiz # center layers and slices for plotting
icy = int(ny/2); cy = ymin + icy*ysiz
icz = int(nz/2); cz = zmin + icz*zsiz

4.5. Kriging para variables continuas - Ley de Cu - UG 4

Aplicamos los conceptos desarrollados anteriormente y el modelo de variograma obtenidos para la UG bajo análisis.

Este capítulo presenta los programas de kriging de GSLIB. Estos programas permiten estimar un volumen ya sea usando kriging simple (SK), kriging ordinario (OK) o kriging con tendencias (KT).

El programa kt3d ofrece un programa de kriging 3D bastante avanzado para puntos o bloques mediante kriging simple (SK), kriging ordinario (OK) o kriging con un modelo de tendencia polinómica (no sera revisado en este curso).

Los parámetros requeridos por el programa principal kt3d se listan a continuación:

• datafl: archivo de entrada con los datos en formato GSLIB.

• DH, X, Y, Z, icolvr, icolsec: columnas para la columna que contiene los codigos DH (ya los borramos y no seran necesarios - valor 0) coordenadas x, y, z, la variable a estimar y la variable de deriva externa (o media no estacionaria - no utilizada en este curso).

• tmin y tmax: triming limits o valores fuera de este rango (menores estrictamente que tmin o mayores o iguales que tmax) se ignoran.: valores fuera de este rango (menores estrictamente que tmin o mayores o iguales que tmax) se ignoran..

• option: 0 para kriging sobre una rejilla de puntos o bloques, 1 para valida validación cruzada con los datos en datafl, 2 para jackknife con los datos del archivo siguiente.

• jackfl: archivo con las ubicaciones para realizar estimaciones (opción jackknife).

• icoljx, icoljy, icoljz, icoljvr, icoljsec: columnas para x, y, z, la variable y la variable secundaria en jackfl.kfl.

• idbg: nivel entero de depuración entre 0 y 3. A mayor nivel, más información entregada. Los niveles usuales son 0 y 1; 2 y 3 imprimen todas las matrices de kriging y datos usados para cada punto/bloque (no recomendable con grillas grandes).

• dbgfl: archivo donde se escribe la salida de depuración.

• outfl: archivo de salida de la rejilla. Contiene la estimación y la varianza de kriging para cada punto/bloque en la rejilla, con orden rápido en x, luego y y finalmente z. Los puntos no estimados se marcan con un número negativo grande (-999.).

• nx, xmn, xsiz: definición de la grilla en la dirección x.

• ny, ymn, ysiz: definición de la grilla en la dirección y.

• nz, zmn, zsiz: definición de la grilla en la dirección z.

• nxdis, nydis, nzdis: número de puntos de discretización para un bloque. Si nxdis, nydis y nzdis son 1, se realiza kriging puntual.

• ndmin y ndmax: número mínimo y máximo de datos a usar para kriging de un bloque.

• noct: máximo a retener por octante (no se usa búsqueda por octantes si noct = 0).i noct = 0).

• radius_hmax, radius_hmin, radius_vert: radios de búsqueda en la dirección horizontal máxima, horizontal mínima y vertical.

• sang1, sang2, sang3: ángulos que describen la orientación de la elipsoide de búsqueda (anisotropía).

• ikrige y skmean: si ikrige = 0 se realiza kriging simple estacionario con media skmean; si ikrige = 1 se realiza kriging ordinario y el valor de skmean no es considerado.

• nst y c0: número de estructuras del variograma y constante nugget. El efecto pepita no se cuenta como estructura.

Para cada una de las nst estructuras anidadas se debe definir el tipo (it), la contribución (cc), los ángulos (azi, dip) y los alcances (a_hmax, a_hmin, a_vert).

Vemos como queda nuestro archivo de parámetros para la UG 4.

parfile = '''                  Parameters for KT3D
                  *******************

START OF PARAMETERS:
./01_data/data_UG4.dat              -file with data
0  1  2  3  4  0                 -   columns for DH,X,Y,Z,var,sec var
-1.0e21   1.0e21                 -   trimming limits
0                                -option: 0=grid, 1=cross, 2=jackknife
xvk.dat                          -file with jackknife data
1   2   0    3    0              -   columns for X,Y,Z,vr and sec var
0                                -debugging level: 0,1,2,3
./03_resultados/kt3d-bl_UG4.dbg                   -file for debugging output
./03_resultados/kt3d-bl_UG4.out                   -file for kriged output
128   471116.24    15                       -nx,xmn,xsiz - AR: pegamos los valores de la grilla definida en la sección anterior
97   6924721.57    15                       -ny,ymn,ysiz
100   2950.02    15                     -nz,zmn,zsiz
4    4      2                    -x,y and z block discretization
12    24                           -min, max data for kriging - AR : uso de 24 datos para el krigeo
3                                -max per octant (0-> not used) - AR : número de muestras por octante
200.0  100.0  90.0                -maximum search radii - AR : radios de búsqueda de acuerdo a alcances de variograma
 45   60   0.0                 -angles for search ellipsoid - AR : ángulos de búsqueda de acuerdo a azimuts y dips principales del variograma
1     0.37                      -0=SK,1=OK,2=non-st SK,3=exdrift - AR : uso de kriging ordinario, segundo valor no es relevante en este caso
0 0 0 0 0 0 0 0 0                -drift: x,y,z,xx,yy,zz,xy,xz,zy
0                                -0, variable; 1, estimate trend
extdrift.dat                     -gridded file with drift/mean
4                                -  column number in gridded file
3 0.003  -nst, nugget effect  - AR: variograma de CU para UG 4  modelado en sección anterior
2 0.022 45 60  0.0    -it,cc,ang1,ang2,ang3
         200 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 9999 90  -a_hmax, a_hmin, a_vert
'''

with open("02_parametros/kt3d_bl_UG4.par", "w") as f:
    f.write(parfile)

Ahora ya podemos ejecutar kt3d de GSLIB. Puede tardar varios minutos dependiendo del tamaño de la grilla definida.

# Ejecutar kt3d de GSLIB
result = subprocess.run(["./Gslib90/kt3d", "./02_parametros/kt3d_bl_UG4.par"], capture_output=True, text=True, check=True)

Se procede a leer los resultados de kriging de indicadores generados por kt3d para la primera UG.

estCU_UG4 = GSLIB.GSLIB2ndarray_3D("./03_resultados/kt3d-bl_UG4.out ", 0, 1, nx, ny, nz)[0][0]

Podemos visualizar distintos cortes de este archivo usando la siguiente función:

idz = 50; idy = 49; idx = 64   # índices para cortes en z, y, x respectivamente. Valores van entre 0 y nz-1, ny-1, nx-1

plt.subplots()
plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
GSLIB.pixelplt_st(estCU_UG4[idz,:,:],xmin,xmax,ymin,ymax,15,0.0,1.0,'Estimacion de Cu - UG 4 - Vista en planta','X(m)','Y(Mm)','Cu [%]',cmap)

plt.subplots()
plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
GSLIB.pixelplt_st(estCU_UG4[:,idy,:],ymin,ymax,zmin,zmax,15,0.0,1.0,'Estimacion de Cu - UG 4 - Vista Este','Y(Mm)','Z(m)','Cu [%]',cmap)

plt.subplots()
plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
GSLIB.pixelplt_st(estCU_UG4[:,:,idx],ymin,ymax,zmin,zmax,15,0.0,1.0,'Estimacion de Cu - UG 4 - Vista Norte','X(m)','Z(m)','Cu [%]',cmap)

<matplotlib.image.AxesImage at 0x125d47e1590>

4.5.1. Visualizando nuestro modelo de Leyes

continuous = 'cu_pct'

x_coords = np.linspace(xmin + xsiz/2, xmax - xsiz/2, nx)
y_coords = np.linspace(ymin + ysiz/2, ymax - ysiz/2, ny)
z_coords = np.linspace(zmin + zsiz/2, zmax - zsiz/2, nz)

import plotly.io as pio
# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode
import plotly.graph_objects as go
import plotly.express as px
import plotly.offline as py

fig = go.Figure(data=[go.Scatter3d(
    x=df['X'],
    y=df['Y'],
    z=df['Z'],
    marker=dict(color=df[continuous],
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=df[continuous].quantile(0.95),
                size=2.0, # Esta opcion controla el tamaño de los datos
                colorbar=dict(
                    title='Cu [%]',
                    thickness=20,
                    # Add a border to the colorbar
                    outlinecolor='black',  # Color of the border
                    outlinewidth=2,       # Width of the border
                    bordercolor='white',   # Background border color (if applicable)
                    borderwidth=1         # Background border width
                )
                ),
    mode='markers',
    opacity=1
)])

# otras opciones para controlar el color de los elementos de la figura, borrar hovermode para desactivar etiquetas al pasar el mouse
fig.update_layout(    title_text='Estimación de Cu - UG 4',
scene = dict(
                    xaxis = dict(
                        title='Este',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    yaxis = dict(
                        title='Norte',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white"),
                    zaxis = dict(
                        title='Elevación',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    ),
                    font_family="Times New Roman",
                    font_color="black", hovermode=False
                  )

# ocultar la leyenda/autonombre "trace 0" junto al colorbar
fig.update_layout(showlegend=False)

# definir cortes en x, y, z

icz = 60  # índice de corte en z (puede ajustarse), va entre 0 y nz-1

# Agregar cortes del modelo de UGs con la misma escala de colores de plotly
x_slice = np.flip(np.flipud(estCU_UG4[:, :, icx]), axis=1)
y_slice = np.flip(estCU_UG4[:, icy, :], axis=0)
z_slice = np.flip(estCU_UG4[nz-icz, :, :], axis=0)

# X-slice (plano constante x)
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Y-slice (plano constante y)
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None,
)

# Z-slice (plano constante z)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=df[continuous].quantile(0.95),
    showscale=False,
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Añadir superficies al figure (manteniendo la leyenda y colores consistentes)
fig.add_trace(surf_x)
fig.add_trace(surf_y)
fig.add_trace(surf_z)

fig.show()

4.6. Categorizacion de Recursos - UG 4

Haremos nuestra categorizacion de recursos para la UG bajo analisis de acuerdo a la varianza de kriging, el cual es un resultado de utilizar kt3d.

## Varianza de la estimación viene en la segunda columna (1)           ↓        del archivo de salida de kt3d
varCU_UG4 = GSLIB.GSLIB2ndarray_3D("./03_resultados/kt3d-bl_UG4.out ", 1, 1, nx, ny, nz)[0][0]

# pio.renderers.default = "colab"
pio.renderers.default = "notebook" # usar esta linea en jupyter notebook o vscode

continuous = 'cu_pct'

var_block_model = varCU_UG4 # setiamos el modelo de varianzas

# encontrar el percentil 99 para definir el rango de colores
max_var = np.nanpercentile(var_block_model, 99.9)

x_coords = np.linspace(xmin + xsiz/2, xmax - xsiz/2, nx)
y_coords = np.linspace(ymin + ysiz/2, ymax - ysiz/2, ny)
z_coords = np.linspace(zmin + zsiz/2, zmax - zsiz/2, nz)


fig = go.Figure(data=[go.Scatter3d(
    x=df['X'],
    y=df['Y'],
    z=df['Z'],
    marker=dict(color=df[continuous]*0.0, # usar un color neutro para los puntos // máxima confianza
                colorscale=px.colors.sequential.Rainbow[1:],
                cmin=0.0,
                cmax=max_var,
                size=2.0, # Esta opcion controla el tamaño de los datos
                colorbar=dict(
                    title='Varianza',
                    thickness=20,
                    # Add a border to the colorbar
                    outlinecolor='black',  # Color of the border
                    outlinewidth=2,       # Width of the border
                    bordercolor='white',   # Background border color (if applicable)
                    borderwidth=1         # Background border width
                )
                ),
    mode='markers',
    opacity=1
)])

# otras opciones para controlar el color de los elementos de la figura, borrar hovermode para desactivar etiquetas al pasar el mouse
fig.update_layout(    title_text='Varianza de estimación de Cu - UG 4',
scene = dict(
                    xaxis = dict(
                        title='Este',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    yaxis = dict(
                        title='Norte',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white"),
                    zaxis = dict(
                        title='Elevación',
                        backgroundcolor="white",
                        gridcolor="gray",
                        showbackground=True,
                        zerolinecolor="white",),
                    ),
                    font_family="Times New Roman",
                    font_color="black", hovermode=False
                  )

# ocultar la leyenda/autonombre "trace 0" junto al colorbar
fig.update_layout(showlegend=False)

# definir cortes en x, y, z

icz = 60  # índice de corte en z (puede ajustarse), va entre 0 y nz-1

# Agregar cortes del modelo de UGs con la misma escala de colores de plotly
x_slice = np.flip(np.flipud(var_block_model[:, :, icx]), axis=1)
y_slice = np.flip(var_block_model[:, icy, :], axis=0)
z_slice = np.flip(var_block_model[nz-icz, :, :], axis=0)

# X-slice (plano constante x)
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=max_var,
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Y-slice (plano constante y)
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=max_var,
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None,
)

# Z-slice (plano constante z)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=px.colors.sequential.Rainbow[1:],
    cmin=0.0,
    cmax=max_var,
    showscale=False,
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Añadir superficies al figure (manteniendo la leyenda y colores consistentes)
fig.add_trace(surf_x)
fig.add_trace(surf_y)
fig.add_trace(surf_z)

fig.show()

4.6.1. Descargar mallas de sondaje tipo

from urllib.request import urlretrieve
## Descargar archivos a la carpeta 01_data

urls = [
    "https://raw.githubusercontent.com/arqlm/arqlm.github.io/main/courses/MIN431/_static/Sondajes_mi.dat",
    "https://raw.githubusercontent.com/arqlm/arqlm.github.io/main/courses/MIN431/_static/Sondajes_ii.dat",
]

for url in urls:
    filename = url.split("/")[-1]
    try:
        urlretrieve(url, f"./01_data/{filename}")
        print(f"Descargado: {filename}")
    except Exception as e:
        print(f"Error al descargar {filename}: {e}")

Descargado: Sondajes_mi.dat
Descargado: Sondajes_ii.dat

parfile = '''                  Parameters for KT3D
                  *******************

START OF PARAMETERS:
./01_data/Sondajes_mi.dat              -file with data
0  1  2  3  4  0                 -   columns for DH,X,Y,Z,var,sec var
0   1.0e21                 -   trimming limits
0                                -option: 0=grid, 1=cross, 2=jackknife - AR : usamos la opción 1 de cross-validation
xvk.dat                          -file with jackknife data - AR : no se usa en este caso
1   2   0    3    0              -   columns for X,Y,Z,vr and sec var - AR : no se usa en este caso
0                                -debugging level: 0,1,2,3
./03_resultados/kt3d_limite_mi_UG4.dbg                   -file for debugging output
./03_resultados/kt3d_limite_mi_UG4.out                   -file for kriged output
1   0    15                 -nx,xmn,xsiz - AR : no se usa en este caso debido a cross-validation
1   0    15                 -ny,ymn,ysiz - AR : no se usa en este caso debido a cross-validation
1   0    15                     -nz,zmn,zsiz - AR : no se usa en este caso debido a cross-validation
4    4      2                    -x,y and z block discretization
1    32                           -min, max data for kriging - AR : uso de 8 datos para el krigeo
0                                -max per octant (0-> not used) - AR : 0 si no se desea usar octantes
200.0  200.0  200.0                -maximum search radii - AR : radios de búsqueda de acuerdo a alcances de variograma
 45   60   0.0                 -angles for search ellipsoid - AR : ángulos de búsqueda de acuerdo a azimuts y dips principales del variograma
1     0.37                      -0=SK,1=OK,2=non-st SK,3=exdrift - AR : uso de kriging ordinario, segundo valor no es relevante en este caso
0 0 0 0 0 0 0 0 0                -drift: x,y,z,xx,yy,zz,xy,xz,zy
0                                -0, variable; 1, estimate trend
extdrift.dat                     -gridded file with drift/mean
4                                -  column number in gridded file
3 0.003  -nst, nugget effect  - AR: variograma de CU para UG 4  modelado en sección anterior
2 0.022 45 60  0.0    -it,cc,ang1,ang2,ang3
         200 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 9999 90  -a_hmax, a_hmin, a_vert
'''

with open("02_parametros/kt3d_limite_mi_UG4.par", "w") as f:
    f.write(parfile)

Ahora ya podemos ejecutar kt3d de GSLIB. Puede tardar varios minutos dependiendo del tamaño de la grilla definida.

# Ejecutar kt3d de GSLIB
result = subprocess.run(["./Gslib90/kt3d", "./02_parametros/kt3d_limite_mi_UG4.par"], capture_output=True, text=True, check=True)

## Leemos la segunda columna del archivo de salida de kt3d
var_lim_mi_UG4 = []
with open("./03_resultados/kt3d_limite_mi_UG4.out", "r") as f:
    for i, line in enumerate(f):
        if i >= 4:  # Start from the fifth line
            columns = line.split()
            if len(columns) > 1:  # Check if there is a second column
                var_lim_mi_UG4.append(float(columns[1]))

var_lim_mi_UG4

[0.01997562]

parfile = '''                  Parameters for KT3D
                  *******************

START OF PARAMETERS:
./01_data/Sondajes_ii.dat              -file with data
0  1  2  3  4  0                 -   columns for DH,X,Y,Z,var,sec var
0   1.0e21                 -   trimming limits
0                                -option: 0=grid, 1=cross, 2=jackknife - AR : usamos la opción 1 de cross-validation
xvk.dat                          -file with jackknife data - AR : no se usa en este caso
1   2   0    3    0              -   columns for X,Y,Z,vr and sec var - AR : no se usa en este caso
0                                -debugging level: 0,1,2,3
./03_resultados/kt3d_limite_ii_UG4.dbg                   -file for debugging output
./03_resultados/kt3d_limite_ii_UG4.out                   -file for kriged output
1   0    15                 -nx,xmn,xsiz - AR : no se usa en este caso debido a cross-validation
1   0    15                 -ny,ymn,ysiz - AR : no se usa en este caso debido a cross-validation
1   0    15                     -nz,zmn,zsiz - AR : no se usa en este caso debido a cross-validation
4    4      2                    -x,y and z block discretization
1    32                           -min, max data for kriging - AR : uso de 8 datos para el krigeo
0                                -max per octant (0-> not used) - AR : 0 si no se desea usar octantes
200.0  200.0  200.0                -maximum search radii - AR : radios de búsqueda de acuerdo a alcances de variograma
 45   60   0.0                 -angles for search ellipsoid - AR : ángulos de búsqueda de acuerdo a azimuts y dips principales del variograma
1     0.37                      -0=SK,1=OK,2=non-st SK,3=exdrift - AR : uso de kriging ordinario, segundo valor no es relevante en este caso
0 0 0 0 0 0 0 0 0                -drift: x,y,z,xx,yy,zz,xy,xz,zy
0                                -0, variable; 1, estimate trend
extdrift.dat                     -gridded file with drift/mean
4                                -  column number in gridded file
3 0.003  -nst, nugget effect  - AR: variograma de CU para UG 4  modelado en sección anterior
2 0.022 45 60  0.0    -it,cc,ang1,ang2,ang3
         200 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 100 90  -a_hmax, a_hmin, a_vert
2 0.005 45 60  0.0    -it,cc,ang1,ang2,ang3
         9999 9999 90  -a_hmax, a_hmin, a_vert
'''

with open("02_parametros/kt3d_limite_ii_UG4.par", "w") as f:
    f.write(parfile)

Ahora ya podemos ejecutar kt3d de GSLIB. Puede tardar varios minutos dependiendo del tamaño de la grilla definida.

# Ejecutar kt3d de GSLIB
result = subprocess.run(["./Gslib90/kt3d", "./02_parametros/kt3d_limite_ii_UG4.par"], capture_output=True, text=True, check=True)

## Leemos la segunda columna del archivo de salida de kt3d
var_lim_ii_UG4 = []
with open("./03_resultados/kt3d_limite_ii_UG4.out", "r") as f:
    for i, line in enumerate(f):
        if i >= 4:  # Start from the fifth line
            columns = line.split()
            if len(columns) > 1:  # Check if there is a second column
                var_lim_ii_UG4.append(float(columns[1]))

var_lim_ii_UG4

[0.02679717]

## Transformamos la varianza del modelo de bloques varCU_UG4 en una variable que indique la categoría del bloque usando numpy

categoria_UG4 = np.empty(varCU_UG4.flatten().shape, dtype=object)
categoria_UG4[:] = "Sin Categoría"  # Valor por defecto para bloques sin datos
categoria_UG4[(varCU_UG4.flatten() < var_lim_mi_UG4) & (varCU_UG4.flatten() >= 0)] = "Medido"  # Categoría 1: Varianza baja
categoria_UG4[varCU_UG4.flatten() >= var_lim_ii_UG4] = "Inferido"  # Categoría 3: Varianza alta
categoria_UG4[(varCU_UG4.flatten() < var_lim_ii_UG4) & (varCU_UG4.flatten() >= var_lim_mi_UG4)] = "Indicado"  # Categoría 2: Varianza media
categoria_UG4 = np.array(categoria_UG4).reshape(varCU_UG4.shape)

fig, ax = plt.subplots(figsize=(5,5))
ax.clear()

plt.rcParams['figure.dpi'] = 250 # Setear resolución de la figura
plt.rcParams["font.family"] = "Times New Roman"

# Proporciones de cada categoría en nuestra base de datos ordenados por Inferido, Indicado, Medido
color = ['#d73027', '#fee08b', '#1a9850']   # Colores para Inferido, Indicado, Medido

proportions = pd.Series(categoria_UG4.flatten()).value_counts(normalize=True).sort_index()
labels = proportions.index.tolist()

# barras ordenadas por Inferido, Indicado, Medido
ordered_labels = ['Inferido', 'Indicado', 'Medido']
proportions = proportions.reindex(ordered_labels)

bars = ax.bar(proportions.index, proportions.values, color=color)

ax.set_title('Proporción de cada categoría en la base de datos')
ax.set_xlabel('Categorías')
ax.set_ylabel('Proporción')
ax.yaxis.set_major_formatter(mticker.PercentFormatter(1.0))

for bar, p in zip(bars, proportions.values):
    ax.annotate(f"{p*100:.1f}%", xy=(bar.get_x() + bar.get_width() / 2, p),
                xytext=(0, 3), textcoords="offset points", ha='center', va='bottom')

fig.tight_layout()
plt.show()

# Ensure renderer is appropriate for your environment (notebook / vscode)
pio.renderers.default = 'notebook'

# --- Input assumptions (these must exist in the environment) ---
# df: pandas DataFrame with columns 'X','Y','Z'
# categoria_UG4: numpy array with shape (nz, ny, nx) containing strings: 'Inferido','Indicado','Medido'
# x_coords, y_coords, z_coords: 1D coordinate arrays for grid cell centers (lengths nx, ny, nz)
# icx, icy, icz: integer indices for the central slices to display
# nz, ny, nx: grid dimensions

# Map categories to numeric codes
cat_to_num = {'Inferido': 0, 'Indicado': 1, 'Medido': 2}
cat_num = np.vectorize(lambda x: cat_to_num.get(x, np.nan))(categoria_UG4)
cat_num = cat_num.astype(float)  # allow NaNs

# Create slices (same convention used in the notebook)
x_slice = np.flip(np.flipud(cat_num[:, :, icx]), axis=1)   # X-slice (nz, ny)
y_slice = np.flip(cat_num[:, icy, :], axis=0)             # Y-slice (nz, nx)
z_slice = np.flip(cat_num[nz-icz, :, :], axis=0)          # Z-slice (ny, nx)

# Discrete colors: Inferido (red), Indicado (yellow), Medido (green)
# We use normalized stops [0.0, 0.5, 1.0] so numeric codes 0,1,2 map exactly to the three colors
colorscale = [
    [0.0, '#d73027'],   # Inferido -> 0
    [0.5, "#fff700"],   # Indicado -> 1 (normalized to 0.5)
    [1.0, '#1a9850']    # Medido -> 2 (normalized to 1.0)
]

# Build surface for X-slice
YY, ZZ = np.meshgrid(y_coords, z_coords)             # shapes (nz, ny)
XX = np.full_like(YY, x_coords[icx])
surf_x = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=x_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=False,
    opacity=1,
    name=f'corte X={x_coords[icx]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Build surface for Y-slice
XX, ZZ = np.meshgrid(x_coords, z_coords)             # shapes (nz, nx)
YY = np.full_like(XX, y_coords[icy])
surf_y = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=y_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=False,
    opacity=1,
    name=f'corte Y={y_coords[icy]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Build surface for Z-slice (we enable a colorbar here with categorical ticks)
XX, YY = np.meshgrid(x_coords, y_coords)             # shapes (ny, nx)
ZZ = np.full_like(XX, z_coords[icz])
surf_z = go.Surface(
    x=XX, y=YY, z=ZZ,
    surfacecolor=z_slice.astype(float),
    colorscale=colorscale,
    cmin=0, cmax=2,
    showscale=True,
    colorbar=dict(
        title='Categoría',
        tickvals=[0, 1, 2],
        ticktext=['Inferido', 'Indicado', 'Medido'],
        thickness=20
    ),
    opacity=1,
    name=f'corte Z={z_coords[icz]:.1f}',
    hoverinfo='skip',
    hovertemplate=None
)

# Add sample points as neutral black markers for reference
scatter_samples = go.Scatter3d(
    x=df['X'], y=df['Y'], z=df['Z'],
    mode='markers',
    marker=dict(color='black', size=2, opacity=0.8),
    name='Muestras'
)

# Compose figure
fig = go.Figure(data=[scatter_samples, surf_x, surf_y, surf_z])
fig.update_layout(
    title='Categoría UG4 (Inferido=rojo, Indicado=amarillo, Medido=verde)',
    scene=dict(
        xaxis=dict(title='Este'),
        yaxis=dict(title='Norte'),
        zaxis=dict(title='Elevación')
    ),
    width=900,
    height=720,
    font=dict(family='Times New Roman', size=12)
)

fig.show()

4.7. Exportar modelo

Una vez verificado que las muestras iniciales coinciden hacen sentido con el modelo de unidades geologicas, podemos exportar este modelo de bloques, el cual sera utilizado nuevamente una vez finalizada la estimacion de leyes dentro de cada unidad para compilar el modelo de leyes finales.

## Exportar modelo de UGs en formato GSLIB
GSLIB.ndarray2GSLIB_3D(modelo_ug.astype(int),'./03_resultados/modelo_ug.dat', 'UG')